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We present a quantum Monte Carlo (QMC) technique for calculating the exact finite-temperature 
properties of Bose-Fermi mixtures. The Bose-Fermi Auxiliary-Field Quantum Monte Carlo (BF- 
AFQMC) algorithm combines two methods, a finite-temperature AFQMC algorithm for bosons 
and a variant of the standard AFQMC algorithm for fermions, into one algorithm for mixtures. 
We demonstrate the accuracy of our method by comparing its results for the Bose-Hubbard and 
Bose-Fermi-Hubbard models against those produced using exact diagonalization for small systems. 
Comparisons are also made with mean-field theory and the worm algorithm for larger systems. As 
is the case with most fermion Hamiltonians, a sign or phase problem is present in BF-AFQMC. We 
discuss the nature of these problems in this framework and describe how they can be controlled 
with well-studied approximations to expand BF-AFQMC's reach. The new algorithm can serve 
as an essential tool for answering many unresolved questions about many-body physics in mixed 
Bose-Fermi systems. 
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I. INTRODUCTION 

Ultracold atomic gases loaded into optical traps offer 
the unique possibility of experimentally simulating many 
of the fundamental models of condensed matter physics 
PJ [2] . These systems are clean, and owing to remarkable 
advances in trapping, cooling, and the manipulation of 
inter- and intraparticle interactions, may be studied with 
an unprecedented level of experimental control. One of 
the field's landmark achievements has been the observa- 
tion of the superfluid-Mott insulator transition in Bose 
gases [3J. Analogous successes with fermions have led 
to the direct observation of such phenomena as Fermi 
pressure and antibunching [U [5] . Focus has now shifted 
to ultracold mixtures of bosons and fermions [BlfTZ] . At 
the most practical level, bosons may be used to sym- 
pathetically cool trapped fermions jT3l [14] . Much more 
tantalizing, however, is the prospect that bosons may be 
able to mediate a BCS superfluid transition in ultracold 
Fermi gases [15H17) . or emulate many-body Hamiltoni- 
ans of mixture systems predicted to exhibit a plethora 
of exotic phases [HI [19] . Equally intriguing is the pos- 
sibility of using newly created "Bose-Fermi molecules" 
with permanent dipole moments as qubits for quantum 
computers or as probes of the permanent electric dipole 
moment of the electron (TUl l2"UH2"2"] . These possibilities 
have galvanized both experimentalists and theorists to 
develop new tools capable of exploring the full range of 
mixture phenomenology. 

From a theoretical standpoint, delineating the exact 
finite-temperature Bose-Fermi phase diagram represents 
a formidable challenge. Mean-field and perturbation the- 
ory calculations suggest that Bose-Fermi mixtures may 
exhibit a wide variety of behaviors, ranging from Bose- 
Fermi "molecule" spin and charge density waves to phase 



segregation [TBI 123 I2"3H2"T] . Nevertheless, these tech- 
niques are approximate by definition, which raises con- 
cerns about the phase diagrams they yield. A reliable de- 
scription of Bose-Fermi mixture phenomenology requires 
an exact framework capable of accurately accounting 
for strong correlation among particles. Accurate results 
can be obtained for small clusters whose limited Hilbert 
spaces are amenable to exact diagonalization (ED), and 
linear chains for which quantum Monte Carlo (QMC) 
techniques free of the sign problem or density matrix 
renormalization group methods may be applied |28H33j . 
Techniques for large systems in two and higher dimen- 
sions, however, are scarce. 

The most promising and flexible technique for mix- 
tures to date uses the framework of Dynamical Mean 
Field Theory (DMFT) [31]. While initial applications of 
DMFT to mixtures paired well-established DMFT meth- 
ods for fermions with approximate treatments of bosons 
[35l437j . the first rigorous Bose-Fermi DMFT algorithm 
has recently been proposed, which weds fermion DMFT 
with a newly-derived DMFT approach for bosons [381 - 
|4"U] . As with all DMFT approaches, this technique is 
only expected to be accurate in the limit of large dimen- 
sionality or coordination number. Indeed, recent Boson- 
DMFT (BDMFT) calculations on the Bose-Hubbard 
model demonstrate that, while DMFT is remarkably ac- 
curate in three dimensions, it is less so in two dimen- 
sions [30] • Furthermore, because DMFT is most useful 
for systems with short-range correlations, inhomogeneous 
phases and long- wavelength collective modes may present 
additional challenges. 

In contrast, QMC techniques offer the promise of being 
exact regardless of system size, dimensionality, and ho- 
mogeneity. QMC techniques differ widely in detail from 
algorithm to algorithm, but all employ stochastic sam- 
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pling to solve the Schrodinger equation at zero tempera- 
ture or determine partition and correlation functions at 
finite temperatures. Because of their accuracy and mod- 
est computational cost, QMC methods such as the world- 
line and worm algorithms have become the techniques 
of choice for boson lattice models [4"TM44| . Auxiliary- 
field and diagrammatic QMC techniques also exist for 
fermions [45-50J. Unlike techniques for bosons, however, 
fermion QMC in two or more dimensions is generally 
plagued by the sign problem, resulting in an exponential 
scaling of computational cost with inverse temperature 
to achieve a fixed accuracy [51] . Developing a widely- 
applicable QMC technique for mixtures thus requires not 
only marrying two considerably different fermion and bo- 
son techniques together, but finding a way to tame the 
sign problem within that combined formalism. 

Widely employed in condensed matter and nuclear 
physics, the Auxiliary-Field Quantum Monte Carlo 
(AFQMC) method [12 IS2 |53] is a field theoreti- 
cal method where many-body propagators resulting 
from two-body interactions are transformed into many- 
dimensional integrals over one-body propagators using 
the Hubbard-Stratonovich Transformation [SSI [57]. The 
resulting integrals are then computed using Monte Carlo 
sampling. In recent years, AFQMC has predominantly 
been used to study the equilibrium properties of the 
Hubbard model both at finite-temperature and in the 
ground-state. Like all fermion QMC techniques, conven- 
tional AFQMC suffers from the sign problem in most pa- 
rameter regimes. However, an alternative formulation, 
in which walkers are pruned using population control 
techniques as they sample AFs in imaginary time, has 
allowed a general, efficient approach to treat both lo- 
cal and extended interactions. This framework allows 
the constrained-path and phaseless approximations to be 
easily incorporated to control the sign and phase prob- 
lems [4"?l I58H6T)] . In recent years, these approximations 
have been tested on a variety of systems including the 
Hubbard model [13 [SHI HO] and the electronic structure 
of solids and molecules (ST] and has been shown to 
yield accurate energies and correlation functions. Thus, 
Constrained-Path AFQMC (CPMC) is well-equipped to 
explore phases beyond the scope of other fermion QMC 
methods. The formalism of AFQMC has also previ- 
ously been generalized to treat bosons in the ground state 
[63 [64]. This suggests that AFQMC would be perfectly 
suited for studying mixtures via a combination of bosonic 
and fermionic Monte Carlo techniques if the formalism 
could be further expanded to treat bosons at finite tem- 
peratures. 

In this work, we present an exact QMC methodology 
that can be used to determine the thermodynamic prop- 
erties of Bose-Fermi mixtures in any dimension over a 
wide range of parameters. Our method, Bose-Fermi Aux- 
iliary Field Quantum Monte Carlo (BF- AFQMC), gener- 
alizes finite-temperature AFQMC for fermions to bosons 



and Bose-Fermi mixtures. By casting the bosonic portion 
of the problem in terms of auxiliary fields, we can extend 
determinantal QMC techniques to bosons and sample 
the boson partition function by sampling determinants 
just as one would for fermions. We arrive at an exact 
technique for mixtures by combining our approach for 
bosons with previous AFQMC techniques for fermions. 
We then discuss how the constrained path and phase- 
less approximations can be imposed to remove the sign 
and phase problems in our method. As a benchmark, we 
compare our algorithm's results for Bose-Hubbard and 
spin-polarized Bose-Fermi-Hubbard clusters to those ob- 
tained using ED. We also contrast our results with those 
from mean-field theory (MFT) and the worm algorithm. 

Our paper is organized as follows: In Section II, we 
begin by reviewing the AFQMC formalism for fermions 
as background for our new algorithm. We then proceed 
to present the underlying formalism for our new boson 
and Bose-Fermi algorithms in Section III, including im- 
portance sampling schemes. We also outline the imple- 
mentation of the constrained path and phaseless approx- 
imations, which can respectively control the sign and 
phase problems. In Section IV, we compare our algo- 
rithm's results for the Bose-Hubbard and spin-polarized 
Bose-Fermi Hubbard models against those produced us- 
ing alternative methods in an effort to demonstrate the 
accuracy of our technique. We finally conclude in Sec- 
tion V, leaving the derivation of the expression relating 
the boson partition function to a determinant and other 
details to the appendices. 

II. PRELIMINARIES 

A. GENERIC MIXTURE HAMILTONIAN AND 
DEFINITIONS 

To facilitate the subsequent discussion, we use the fol- 
lowing form of the Bose-Fermi-Hubbard Hamiltonian as 
a concrete example 

H bf =K b + K f + V b + V f + V c , (1) 
where K b contains all one-body boson terms 

Kb = -tb E ffii + H - c ) + E ( 2 ) 

(ij) < 

Kf contains all one-body fermion terms 

K S = -*f E (4/^ + ^-)+E e {>^ ( 3 ) 

(ij),a i,cr 

V b contains two-body boson terms, 
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Vf contains two-body fermion terms 

V f = UfY^m^m^ 

i 



(5) 



and V c represents the Bose-Fermi coupling term 

V c = C'^hirh i . (6) 



the short-time propagators 
yields 



At second order this 



e -Ar(K f +V f ) _ e -(l/2)ATK f e -ArV f e -(l/2)Ark f + q, At 3^ 

(9)' 

which becomes exact in the limit At — > 0. Each short- 
time propagator is thus a product of two one-body prop- 
agators and one two-body propagator. In our Hamilto- 



In the above, b\ , hi denote the boson creation and anni- 
hilation operators and /? , fi„ the fermion creation and 
annihilation operators with spin a (=f or i) a t site i. 
We define the boson density at site i as hi = b\bi and 
the fermion densities as = f} a fi a . The total fermion 
density at each site is denoted by = m^-f -I- m^. % 
and tf represent the respective boson and fermion hop- 
ping parameters. Ub is the two-body boson-boson poten- 
tial, Uf is the two-body fermion-fermion potential, and 
C is the Bose-Fermi coupling. e\ and e{ represent coeffi- 
cients of one-body terms that may include contributions 
from chemical potentials, external traps, or disorder. De- 
pending upon the values of the various parameters, this 
Hamiltonian can exhibit the full range of Bose-Fermi phe- 
nomenology. More general Hamiltonians may be handled 
by the approach outlined below. 



B. FINITE-TEMPERATURE AFQMC FOR 
FERMIONS 



The finite-temperature AFQMC method for fermions 
calculates the thermodynamic properties of a system of 
particles with two-body interactions by reexpressing two- 
body propagators as integrals over one-body propagators 
and a set of auxiliary fields. Here, we review the ba- 
sic formalism to acquaint the reader with previous work 
relevant to the following discussion [65]. In general, the 
finite-temperature expectation value of an observable, O, 
may be written as 



(6) = 



Tr { 6e~ pfl 



Tr 



(7) 



where H is the Hamiltonian of the system and (3 — 
1/kgT. One may rewrite the partition function, Z, in 
terms of a product of I short-time propagators 



Z = Tr 



(e-^) = Tr ( 



e ^H e -ArH_ 



-AtH 



(8) 



Here, At = [3/1 is the timeslice in imaginary time. 
For simplicity, consider the fermion Hamiltonian Hf = 
Kf + Vf using the definitions from Part A. One may 
next perform a Trotter-Suzuki factorization on each of 



Uf 
2 



E 



(10) 



U£ V 
2 ^ 



[m 



it 



hi 



This form allows the two-body propagators to be reex- 
pressed in terms of an integral over a product of one- 
body propagators and a set of auxiliary fields using the 
Hubbard-Stratonovich (HS) Transformation [5T>1 157j : 



,(l/2)Ar« 2 



1 



/2tt 



d(f>e 



"(1/2) 



o 0\/Arii 



(ii) 



where </> is an auxiliary-field (AF). Note that, while there 
are discrete versions of the HS transformation for the 
form of Vf in our Hamiltonian, we have outlined a con- 
tinuous version which formally resembles the transforma- 
tion we will use for Vb and V c . 

This expression for the short-time propagator may be 
further simplified by viewing the collection of fields at 
each timeslice as a vector of fields, <j> = {fa, fa, ...,fav}, 
where N is the number of lattice sites, and the normalized 
Gaussians at each site as probabilities, p(fa). Collecting 
all one-body operators into Bf(fa), we arrive at [6"S"] : 



e -\ArK s e -ArV s e -\ArK s = 



d$ P (faBf(fa, (12) 



where 



B f (fa 



-\AtK 



(rhif-rhn) 



e~i ATK f. 



(13) 

and the one-body term in Eq. [10] can be absorbed by 
replacing e{ a with (e{ a + Uf/2) in Eq. |i] 

Substituting Eq. [Tjjinto the expression for Z in Eq. [HJ 
one arrives at the central AFQMC equation 



d$ P ($)Tr[Bf(fa)...Bf(fa 



(14) 



where $ denotes the full collection of auxiliary-fields at 
each timeslice and site and p($) is the corresponding 
probability of selecting those fields. 

The partition function may therefore be viewed as an 
integral over all fields of the Gaussian probability of se- 
lecting a set of fields multiplied by the trace of single- 
body operators evaluated as a function of the fields. The 
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set of fields at each timeslice and site constitutes a path 
in AF space. Thus, in AFQMC, one calculates the multi- 
dimensional partition function by stochastically sampling 
a set of paths in AF space and evaluating the weighted 
average of the trace along those paths. 

It turns out that the fermion trace over one-body prop- 
agators can be evaluated analytically and expressed as a 
determinant 1691 



Tr f ( 



B 



fVPl 



)-Bf{4>x) 



Det 



I + B f ( 



■Bt 



(15) 

If the size of the single-particle basis (in this case the 
number of lattice sites) is N, Bf(cf>k) is an NxN matrix 
of the propagator Bf(<fik) expressed in that basis, and / is 
the corresponding unit matrix. Inserting this expression 
into that for the partition function, one arrives at 



Zf = 



d<$>p{<$>)Det 



B 



f{<Pl 



)-B f ( 



(16) 



In a similar vein, tracing over fermionic operators yields 
the fermion Green's function: 



Tr 



G 



f (f i f}B f ($ l )...B f ($ 1 )) 

Tr f (Bj ($,).. .Bfifa) 
I 

I + Bj${)...B f {${) 



(17) 
(18) 



where the subscripts on the right denote the (i, j)th el- 
ement of the matrix. Most observables of interest may 
be easily expressed in terms of the single-particle Green's 
function using Wick's theorem [68j . 

With Eqs. [T6| and [T8| in hand, one can evaluate nearly 
any observable by sampling paths according to the parti- 
tion function, calculating the Green's function (and hence 
any related observable) as a function of those paths, and 
weighting the resulting values by the probability of the 
paths sampled. We next present the formalism that al- 
lows one to do the same for bosons. 



III. METHODS 



As we show in Appendix A, the trace over bosons may 
also be expressed as a determinant (which has been noted 
in other contexts before I72H741 ) : 



Tr b [B b (ip l )...B b (ip 1 )) =Det 



J-£ b (^)-5(,(^i) 



allowing the partition function to be expressed as 



d^p($)Det 



I 



I — B b {ipi),.,Bb{ij}i) 



(20) 



.(21) 



Further manipulations yield the boson single-particle 
Green's function 



G 



b _ 



(22) 



Tr (b l b]B b ^ l )...BSi)) 

Tr (B b (^)...B b {ifi)\ 
I 

- B h $ l )...B h $ i ) 



In a boson Auxiliary-Field Quantum Monte Carlo 
(B-AFQMC) algorithm, one can therefore calculate bo- 
son observables by sampling paths according to the bo- 
son partition function in Eq. 21 and evaluating the 



weighted average of observables determined from the bo- 
son Green's function in Eq. |22| There are only two formal 
differences between B-AFQMC and standard fermion 
AFQMC: the minus sign in front of the product of the 
one-body propagators, and the inverse in the determi- 
nant. These differences, however, have a large impact on 
how the B-AFQMC algorithm is implemented compared 
to standard AFQMC. As discussed in detail in Appendix 
B, the new form of the Green's function requires that 
adjustments be made to the way one stabilizes products 
of one-body matrices at low temperatures, while the new 
form of the determinant requires that adjustments be 
made to the way local-updates to the Green's function 
are computed and weights are accumulated as fields are 
selected at each timeslice and site. Except for these ad- 
justments, B-AFQMC maps formally and directly onto 
previous AFQMC algorithms. 



A. FINITE-TEMPERATURE AFQMC FOR 
BOSONS 



B. BOSE-FERMI AFQMC 



Following the same steps outlined for the fermion 
Hamiltonian, Hf, in Section II, one can similarly derive 
an expression relating the boson partition function to in- 
tegrals over one-body boson propagators, B b {il>k), and 
auxiliary fields ip k = {V>ifc, ip2k, 4>Nk}: 



Z b = Tr, 



(e-W») (19) 
' d$p(*)Tr b (bS 1 )...bSi)) ■ 



To combine AFQMC and B-AFQMC into a procedure 
for mixtures, one needs to decouple the Bose-Fermi cou- 
pling term in Eq. [T] This can be done by reexpressing 
Eq. [6] in a form suitable for the HS Transformation: 



V r = 



c 

2 



(23) 



where for brevity we have assumed spin-polarized 
fermions (a =f only). The more general case can be 
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handled similarly by combining the resulting fermion in- 
teraction term with Vf. One may now apply the HS 
Transformation of Eq. [TT]to write each square into linear 
forms as we have shown in Sec. II. B. Note that the result- 
ing h 2 terms can be absorbed into the two-body boson 
term, V b , in Eq. [I| 

An important way to improve the efficiency of BF- 
AFQMC simulations is to subtract any background terms 
prior to the HS Transformation. In both boson and 
fermion ground-state calculations, this was shown to 
greatly reduce the QMC statistical fluctuations and the 
severity of the sign and phase problems [5TJ [53] . For 
example, in Eq. 



23 



one would rewrite (ri. 



v 2 = (v- (v)) 2 + 2v(v) - (v) 2 
= v' 2 + 2{v)v- (v) 2 



(24) 



for each site i, where (v) = {hi + rhi) = {hi) + {rhi), 
with {hi) and {rhi) the average (or desired) boson and 
fermion site densities, e.g., from MFT or exact symmetry 
properties. The HS Transformation is then applied to 

- 2 

v' instead of v 2 , and the one-body and constant terms 
in Eq. 24 can be easily combined with other one-body 



terms in the Hamiltonian and absorbed into the resulting 
one-body propagators, B. 

The background subtraction is intimately connected 
with the mean-field formalism [M]. The idea is to use a 

form of HS Transformation to decouple v' terms which 
are zero in some mean-field framework. That is, setting 
the AF value to zero in the HS decomposition would give 
the corresponding mean- field result. The background 
subtraction is applied to all V b and V c terms; no back- 
ground subtraction is applied to Vt because we have used 



a spin-decomposition (as opposed to charge) in Eq. 10 for 
fermions. The values of {hi) and {rhi) ar e set prior to the 
simulation. It should be emphasized that the formalism 
is exact independent of the choice of mean-field values; 
only the statistical errors are affected. 
The combined partition function is 



Z bf = Tr b 



Tr 



(25) 



After the HS Transformation, the fermion and boson 
propagators are decoupled at each timeslice and site. Be- 
cause all fermion operators commute with all boson oper- 
ators, the propagators may be separated into completely 
independent products of one-body boson and fermion 
propagators. One may then evaluate the traces over these 
products individually to obtain 



Z bf = / d*c?$p(*,$) 

J —OO 

Drt. — — 

I - B b (<ip l )...B b (ip 1 ) 



Del 



(26) 



Because Eq. [T] contains three terms quadratic in the 
boson and fermion densities, three HS Transformations 
must be used at each timeslice and site to reduce these 
terms to one-body operators. The boson and fermion 
Green's functions may analogously be written as above, 
but with one-body matrices that now contain their re- 
spective contributions from the coupling terms. Thus, 
in BF-AFQMC, a generic Bose- Fermi Hamiltonian may 
be simulated by first rewriting all coupling terms such 
that they can be transformed into independent boson and 
fermion propagators. Once the propagators are reparti- 
tioned, the individual boson and fermion Green's func- 
tions may then be evaluated as if there was no coupling 
term, so long as paths are sampled from the full Bose- 
Fermi partition function. 



C. IMPORTANCE SAMPLING 

Determinants are computed using a set of walkers 
whose weights and Green's functions are determined as 
each field is sampled sequentially in imaginary time. 
At the beginning of our simulations, we initialize the 
weights, W($>, ^), of a collection of walkers to 1. Wc 
similarly initialize each walker's Green's function to that 
corresponding to a trial Hamiltonian, such that 



and 



U i3 - 



G 



I-BT...BT 



I + BJ...BJ 



(27) 



(28) 



where B T is a trial one-body matrix at each timeslice. In 
the work that follows, the trial Hamiltonian is typically 
the exact Hamiltonian minus any terms quadratic in the 
density (v 12 terms, after background subtraction). Since 
the chemical potential corresponding to some desired fill- 
ing differs between the trial and exact Hamiltonians, care 
must be taken to determine the appropriate chemical po- 
tential for the trial Hamiltonian before sampling proceeds 
so as to prevent additional statistical fluctuations. 

As each field (or fields, if multiple HS Transforma- 
tions are performed) is selected at site i and timeslice 
k, the weights of the walkers are multiplied by a fac- 
tor, W(<fiik,ipik)- In the absence of importance sampling 
(see below), W{<pik,ipik) is the ratio of the product of 
the newly-updated determinants to the old determinants. 
Let Pf. denote the fermion determinant constructed of 
fields sampled up to the i-th site and fc-th timeslice 



P 



ik 



Del 




■>ik ■ 



■ B f 



(29) 



G 



and P^ k define the corresponding boson determinant 



pb 



Det 



(30) 

where the yet unspecified AF's in the fc-th time slice (for 
sites i through N) can be thought of as having value zero, 
as mentioned in Sec. III.B above. Then, the weight may 
be defined as 



pi pb 

r ik r ik 



P, 



pb 



(31) 



The final product of these factors over all sampled fields 
is proportional to the product of boson and fermion de- 
terminants for the full path that we wish to sample. As 
each field is sampled, the Green's functions are also up- 
dated by replacing the trial one-body matrices with the 
exact one-body matrices based upon the fields. The cor- 
responding Green's function matrix, after sampling field 
i at timeslice fc, would therefore be 



and 



G" 



I 



I - B b T ...B b T B b (^ k ...^ lk )...B b (^, 1 ) 



I + Bj...BjB f ( 



>ik- 



>\k) 



.B 



(32) 



(33) 



f\ 



All trial matrices are replaced until all fields are sampled 
and the Green's functions correspond to those for the 
exact Hamiltonian. After all fields are sampled, aver- 
age observables are computed. The weights and Green's 
functions are then reinitialized to their starting values 
and fields are sampled again until the desired number of 
samples have been collected. 

Of course, if the fields are dra wn randomly according 

will cancel in successive 



31 



to p(<&, 4"), the ratios in Eq. 
steps, and our sampling procedure above will be identi- 
cal to simply sampling entire paths of AFs randomly and 
then calculating the determinants in Eq. [26] as weights 
of the paths. The advantage of the sampling scheme 
above is that it allows importance sampling to be done 
efficiently and, as we discuss in the next subsection, con- 
strained path and phaseless approximations to be easily 
incorporated to control sign and phase problems [2H El] . 

Importance sampling uses an estimated contribution 
based on a trial wave function or density matrix to guide 
the sampling of AFs [?71 [551 • Just as gains in effi- 
ciency may be obtained by subtracting the average den- 
sity from the exact density in each HS-Transformed prop- 
agator, even further gains may be obtained by subtract- 
ing a site-dependent shift, from the auxiliary-field, 
ipi. This shift, called a force bias, effectively modifies the 
probability p(^) for sampling ipi k , to take into account 



the AF paths that have been built up so far, i.e., the prior 
ip values (from ipn to ip(i-i)k)- The shift is added by 
performing a change of variable in the usual HS Trans- 
formation. For example, the boson 2-body term, after 
absorbing the contribution from V c and background sub- 
traction, can be written as 



-&T/2(u b -c){tu-(iu)y 



(34) 



1 



dtpie 



\J2lT j-oo 

^iji e (i>i-j>i)y/-AT(U b -C)(ni-(ni)) 

dil)ip(ipi)W' (ip i ,'ijj i )B(ip i - <4>i) 

where shift- and field-related constants may be regrouped 
into an additional weighting term, W (tpi,ipi), that con- 
tributes to Eq. [31] The one-body operator is now also a 
function of the shift. 

Optimal importance sampling is achieved when the 
shift is chosen such that the fluctuations in the weights 
of the walkers are minimized. At finite-temperatures (see 
the ground-state derivation in Purwanto and Zhang |63|). 
the optimal shift may be shown to be 



ipi 



Tr 


ViB{<${) 


...B$i) 


Tr 







(35) 



where i>i represents the coefficient of the field in the 
HS-Transformed propagator. In the case of Eq. 34 
%i = \] — Ar([/ b — C){hi — (hi)). Shifts may be calcu- 
lated in this way for each HS Transformation. This im- 
portance sampling technique enables us to simulate well 
into the moderate-coupling regime with high-efficiency, 
free of any approximations. 



D. THE CONSTRAINED PATH AND 
PHASELESS APPROXIMATIONS 

As alluded to earlier, a phase problem develops when- 
ever complex propagators produce complex determi- 
nants. When sampled by walkers, these complex deter- 
minants in turn yield complex walker weights. Although 
background subtraction and importance sampling, as dis- 
cussed above, can help reduce statistical fluctuations, the 
phase problem will eventually overwhelm any simulation 
at sufficiently low temperatures or sufficiently large re- 
pulsive interactions. The signature of the sign or phase 
problem is that the weights will populate both positive 
and negative values on the real axis (sign problem) or 
arbitrary phase angles in the complex plane, resulting in 
dramatic cancellation and large fluctuations. The phase 
problem may be avoided with the phaseless approxima- 
tion, an approximation that renders the weights of com- 
plex walkers real via a gauge transformation using a trial 
wave function or density matrix |47l 163) . 
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In the phaseless approximation, one first uses impor- 
tance sampling as described in Section III.C to minimize 
the phase of the weighting factor at each step (timeslice 
and site). Without importance sampling, the weighting 
factor is given by Eq. |31| With importance sampling, it 
becomes 



p/ Jjb 

r ik r ik 



:P?. 



-W {(pik,(f>ik,ipik,ipik) -(36) 



W{<f>ik,ipik) = - 

* {i-l)k* (i-l)fc 

With the optimal choice of force bias, as we discussed 
in Eq. [35] it can be shown that the overall phase ac- 
cumulation is proportional to Ar/m(£Jx,), where El is 
the so-called local energy [471 E3]- In the case of the 
exact trial wave function or density matrix, the imagi- 
nary part of El vanishes. Once the phase is optimally 
reduced, the phaseless approximation omits the overall 
phase. It then projects the random walk to the real axis 
to constrain the overall phase to one gauge choice. In 
the finite-temperature phaseless approximation, we de- 
fine the phase rotation angle AO as 



A6 = Im In 




(37) 



The phase angle may more generally be defined in terms 
of the ratios of both the boson and fermion determinants, 
however we find that the phase problem may typically 
be attributed to boson fluctuations in the Hamiltoni- 
ans studied here. We then multiply the modulus of the 
weighting factor, \W((f>ik,ipik)\ by if |A0| > ir/2 and 
cos(A9) otherwise. This keeps the walker weights real, 
preventing the mass cancellation of weights symptomatic 
of a bad phase problem. 

In addition to the phase problem from bosons, a mix- 
ture simulation may also encounter the sign problem for 
fermions at low temperatures [5T] , which is a special case 
of the phase problem. The phaseless approximation in 
the case of a real HS Transformation and real deter- 
minants reduces to the constrained-path approximation 
[58j . We use this approximation to curb the sign problem 
in this situation. As soon as a walker's fermion determi- 
nant becomes negative, its weight is set to zero. We thus 
sample only those paths such that 



Det 



' l-k 

11*? 



> 



(38) 



for all k from to I. As previously discussed in the liter- 
ature, this prevents corrupted paths whose determinants 
have changed sign from contributing to observable aver- 
ages. 

IV. RESULTS 

In this section, we present illustrative results from our 
Bose-Fermi AFQMC method. Results are compared to 



those obtained from ED, MFT, and the boson worm al- 
gorithm [4TJ|42]. Except where indicated, our B- AFQMC 
and BF- AFQMC calculations were done without impos- 
ing the phaseless or constrained path approximations; 
some were done without importance sampling for bench- 
marking or testing purposes. No optimization was per- 
formed on the choice of the parameters such as the Trot- 
ter step and the intervals with which population control 
[55] or stabilization procedures are applied, except to en- 
sure that the resulting bias is well within statistical er- 
rors. 

ED is a method in which exact expectation values are 
calculated from eigenvalues obtained by diagonalizing the 
system Hamiltonian |66) . In the grand canonical ensem- 
ble, one must determine these eigenvalues for all fermion 
and boson particle numbers. Since a system may in prin- 
ciple be occupied by an infinite number of bosons, an 
exact ED answer would require diagonalizing an infinite 
number of canonical ensemble Hamiltonians. In the re- 
sults that follow, we only include a truncated number of 
bosons sufficient to converge our results to within three 
decimal places. Where the system does not collapse, this 
is sufficient. Near collapse, however, the truncation error 
was visible when compared with the BF- AFQMC results 
and it was necessary to increase the number of bosons 
included in the ED. In our simple implementation, only 
small clusters of up to about five lattice sites could be 
converged to the desired filling with this accuracy. 

For larger systems for which ED fails, we compare to 
MFT. MFT results are expected to be accurate only in 
the weak-coupling regime. Nevertheless, they provide a 
check on our results and demonstrate for which parame- 
ters our exact approach should be particularly valuable. 
In our mean field calculations, we use the general Hamil- 
tonian 

H MF = Kb + K f 

+ -y ^ (2fii(n t ) - hj - (fij) 2 ) 
i 

i 

+ C^(hi(mi) +mi(hi) - (ni){mi)) , (39) 

i 

keeping only the appropriate terms for the given model. 
In these calculations, we self-consistently solve for the 
exact boson and fermion densities at each site until our 
answer is converged to within three decimal places. 

Outside of the weak-coupling regime, we compare our 
results for the Bose-Hubbard Model to those obtained 
from the ALPS Projects' implementation of the worm 
algorithm [67] . The worm algorithm yields exact results 
for bosons for any system size, in any coupling regime 
[HI In all of our worm calculations, we capped the 
number of bosons at each lattice site at a value sufficient 
to achieve convergence in the energies and densities. 
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A. BOSE-HUBBARD MODEL 

We begin by benchmarking our results for the Bosc- 
Hubbard model. The Bose-Hubbard model has long 
been the model of choice for studying condensed He-4 in 
porous media |75j . It has recently been revived to model 
ultracold bosons in optical lattices [3J. The Hamiltonian 
is a special case of Eq. [T] with the fermion constants all 
set to 0. For Ub < in Eq. [4] and sufficiently low tem- 
peratures and high densities, the Bose-Hubbard model 
is expected to exhibit collapse [63l [64] . In the examples 
that follow, we therefore only present results for repul- 
sive Ui,. Our results are equally accurate for Ub < 
before the collapse point, however. Since using Ub > 
results in a phase problem, all of the results that follow 
are averaged over complex phases, without the phaseless 
approximation. The QMC results are thus expected to 
be exact. 

As a first check, we consider a 3 x 1 lattice, with 
tb = 0.01, and (n&) = 1. In Figs, [I] and [2j we com- 
pare our results to those from ED for the energies and 
condensate fractions for varying Ub down to tempera- 
tures T/t w .3. Condensate fractions measure the frac- 
tion of the system lying in the lowest eigenstate of the 
Hamiltonian [331 EH] • As we see in both Figs, [I] and |J 
QMC is exact within error bars well beyond where the 
condensate fraction asymptotes to 1. This suggests that 
our technique can calculate correct expectation values 
from high temperatures corresponding to the Mott insu- 
lating regime to low temperatures corresponding to the 
finite-size version of a superfluid. In Fig. [2j we also plot 
the MFT results for the condensate fractions to illus- 
trate the effects of fluctuations. Only one curve is shown 
for the MFT condensate fractions because they are in- 
dependent of Ub/t. It is evident from this figure that 
MFT yields poor approximations to the true condensate 
fractions even at relatively high temperatures and low 
coupling strengths. Indeed, it only reproduces the exact 
condensate fractions throughout this limited temperature 
range for Ub/t = .5. As illustrated below in Fig. |4j even 
in situations where mean-field condensate fractions are 
nearly exact, energies produced using MFT may be un- 
reliable. This underscores the importance of using exact 
methods where possible. 

The data in Figs. [I] and [2] were calculated without im- 
portance sampling or the phaseless approximation. In 
Fig. [3j we show that we obtain the same results with 
improved statistics using these techniques for Ub/t = .5. 
Using importance sampling and the phaseless approxi- 
mation, our error bars on the number of bosons for the 
same number of samples are at least halved compared 
to those obtained without importance sampling. Error 
bars on other quantities are too small to judge. In previ- 
ous works, importance sampling was observed to greatly 
reduce the error bars in finite-temperature fermion cal- 



culations [SB]. Similarly, in ground state boson calcu- 
lations, an order of magnitude or more improvement in 
efficiency is seen [63] [64] . Our phaseless calculations for 
finite-temperature bosons therefore do not see the dra- 
matic error bar reductions seen in other applications. 
There are several reasons for this. The system size is 
small, such that the variations in the sampled space are 
much reduced compared to larger systems, where the ef- 
fect of importance sampling is expected to increase sig- 
nificantly. The present boson finite-temperature calcu- 
lations are performed in the grand canonical ensemble, 
which could contribute to increased fluctuations. The 
main contribution to the statistical fluctuations in the bo- 
son calculations is likely from the so-called "rogue eigen- 
value" problem which we discuss in the next section. 

For larger lattices, we compare to the worm algorithm. 
Fig. [4] demonstrates that B-AFQMC energies are consis- 
tent with worm energies for 2D systems of varying sizes 
for several Ub- Interestingly, as alluded to above, QMC 
and MFT energies differ dramatically at all but the high- 
est of temperatures. This is even so when the energies 
are normalized to account for the fact that the QMC and 
MFT algorithms require different chemical potentials to 
achieve the same fixed boson number. Fig. [4] may read- 
ily be extended to larger lattices and boson-boson repul- 
sions, but at the price of the increased sampling needed 
to surmount the phase problem. 

B. SPIN-POLARIZED BOSE-FERMI-HUBBARD 
MODEL 

In order to illustrate our Bose-Fermi AFQMC method, 
we similarly apply our technique to the Bose-Fcrmi- 
Hubbard model, the standard model for studying ul- 
tracold mixture phenomenology. As mentioned before, 
here, we limit ourselves to the spin-polarized Hamilto- 
nian, namely Eq. [ljwith fhn — 0. 

As with the Bose-Hubbard model, collapse is antici- 
pated for Ub < and any value of C for densities suffi- 
ciently large that the boson-boson attraction term domi- 
nates the linear coupling term. If Ub — and the boson- 
boson interaction does not dominate, collapse may also 
be observed for a sufficiently large and negative C. The 
phase problem is observed whenever C > or Ub > C. 
We thus again simulate amidst the phase problem so as 
to at once avoid collapse and demonstrate the accuracy 
of our algorithm despite complex phases. 

As our first example, we consider a 2-site Bose-Fermi- 
Hubbard model with varying Ub — C, tb — tf = 0.01, 
and (rib) = ( n ff) = 1- We find that our results for the 
potential energies, kinetic energies, condensate fractions, 
and double occupancies per site [71] agree with ED to 
within small error bars for Ub/t = C/t values up to 13. 
Ub/t = C/t ratios up to 7 are shown in Fig. [5] for the sake 
of clarity. These results demonstrate the correctness of 



9 



our algorithm and implementation, and that exact com- 
putations are feasible for moderate coupling strengths 
amidst an appreciable phase problem. We expect that 
our ability to calculate observables amidst such large 
phase problems will diminish with larger system sizes 
where fewer samples may be taken within a fixed time. 
More sophisticated sampling techniques, better handling 
of the "rogue eigenvalue problem" (see below), and the 
use of the phaseless approximation will drastically im- 
prove the statistical accuracy. 

Lastly, as a check on our mixture algorithm for larger 
systems sizes, we compare to results from MFT in the 
limit of small Ub and C. Our results in Fig. [6] are in 
concurrence with those from MFT for up to 8x8 sys- 
tems (larger sizes are not pictured here). A similar 
comparison, not presented here, was made for the Bosc- 
Hubbard Hamiltonian and yielded analogous results. In 
both cases, MFT results compare well with QMC results 
until the two begin to deviate at lower temperatures, as 
expected. Because there are a limited number of exact 
methods for multidimensional mixtures to which we can 
compare, we reserve further mixture examples and appli- 
cations for a future publication. 



V. DISCUSSION 
A. CHALLENGES 

As the results presented in this work demonstrate, our 
algorithm represents, in principle, an exact method for 
simulating the thermodynamic behavior of an essentially 
arbitrary lattice system composed of interacting bosons 
and fermions. Nevertheless, its performance is still hin- 
dered by several practical challenges. 

One of the more benign challenges relates to the es- 
timation of the correct chemical potentials for desired 
fillings. In order to simulate a mixture with the desired 
fillings in the grand canonical ensemble, one must es- 
timate not only the correct fermion chemical potential, 
but the correct boson chemical potential as well. This 
task is particularly laborious for bosons since their fill- 
ings may change especially rapidly with chemical poten- 
tial. When fillings change more gradually with chemical 
potential, such as in the Mott insulator or normal liquid 
regimes, iterative methods may be employed. Outside of 
such regimes, particularly near or in superfluid phases, 
such methods fail because incorrect or unphysical chemi- 
cal potentials may yield seemingly correct fillings within 
error bars. 

A second challenge to our algorithm is posed by the 
phase problem. As discussed in Section III.D, whenever 
propagators become complex, walker weights and Green's 
functions acquire a complex phase. When this phase 
grows particularly large, controlling statistical fluctua- 
tions becomes a computational challenge. The severity of 



the phase problem depends upon the model and simula- 
tion parameters. For the Bose-Hubbard model, the phase 
problem develops for positive Ub', for the Bose-Fermi- 
Hubbard model it is present whenever C > or Ub > C. 
As with the related sign problem in fermion QMC, the 
severity of the phase problem grows exponentially with 
system size or inverse temperature. This means that for 
large systems and at low temperatures, we need to prop- 
erly impose constraints that systematically bias the re- 
sults. The performance of the constraint in ground state 
calculations should provide a "lower bound" to the qual- 
ity of the approximation in these finite-temperature cal- 
culations. As was previously discussed, importance sam- 
pling can significantly reduce statistical fluctuations, and 
where importance sampling fails, the phaseless approxi- 
mation may be invoked. However, how the approxima- 
tion performs across a phase transition, especially when 
the constraining trial density matrix is poor, remains to 
be studied. 

Perhaps the biggest issue in the present formulation 
relates to the fact that in the grand canonical ensemble 
boson numbers may fluctuate in an unbounded manner. 
In the auxiliary-field formalism, the many-body problem 
is turned into multiple independent-particle problems in 
external fields. By fluctuation of the external fields, the 
target chemical potential may be too high for a partic- 
ular independent-particle path, which would result in a 
condensate with an infinite number of particles. We have 
termed this the "rogue eigenvalue problem." 
As seen in Eq. 



20 



our boson partition function is ex- 
pressed as a determinant of a matrix whose denominator 
may approach or fall below zero. This happens when- 
ever the largest eigenvalue of the product of one-body 
boson matrices approaches or surpasses unity. Although 
it is unphysical for the leading eigenvalue to surpass one - 
and indeed, it never does in our completely deterministic 
mean-field calculations - our walkers may stochastically 
sample such unphysical paths and their related "rogue 
eigenvalues." Walkers whose eigenvalues have surpassed 
one at any point in imaginary time possess corrupted 
paths that develop appreciable phase problems more se- 
vere than those seen in fermion systems and unique to 
simulations of bosons in the grand canonical ensemble. 
This is the leading challenge which impacts the effec- 
tiveness of the algorithm even in the presence of impor- 
tance sampling and the phaseless approximation. In or- 
der to obtain sensible results well into condensed phases 
where eigenvalues may approach one on physical grounds, 
we must therefore prevent walkers from sampling rogue 
paths. One facile method for suppressing rogue paths 
used to produce many of the figures in this paper in- 
volved using larger (uj) values than the mean-field values. 



Instead of setting (vi) in Eq. 24 to the sum of the mean- 
field densities at a given site, we set it to larger values 
that increase the effective chemical potential seen by the 
Green's functions. This reduces the risk of a rogue eigen- 
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value problem at the cost of increased phase fluctuations, 
which can be surmounted by increased averaging. Fur- 
ther details about this approach and more sophisticated 
ones will be presented in an upcoming publication. 



B. CONCLUSIONS 

In this work, we have outlined a new algorithm that 
enables the exact calculation of the thermodynamic prop- 
erties of BF mixtures in multiple dimensions over a wide 
range of parameters. This algorithm enables us to sam- 
ple the boson partition function and calculate boson ex- 
pectation values much as one would sample the fermion 
partition function and calculate fermion expectation val- 
ues using conventional fermion AFQMC. Our method is, 
in principle, exact and we have demonstrated its accu- 
racy by comparing our results to those obtained via ED 
and MFT for the Bose-Hubbard and spin-polarized Bosc- 
Fermi Hubbard models. Approximations need only be 
invoked when stochastic errors stemming from the sign 
and phase problems become uncontrollable. Because our 
algorithm is at once exact and computationally tractable, 
we believe it is uniquely positioned to answer many open 
questions about the Bose-Fermi phase diagram and re- 
cent mixture experiments. Our algorithm is particularly 
well-suited for the study of inhomogenous phases with 
long-range correlations, which cannot be reliably cap- 
tured by mean-field approaches. We leave applications 
of our method to problems of genuine physical interest 
to future publications. 
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APPENDIX A: DERIVATION OF BOSON 
PARTITION AND GREEN'S FUNCTIONS 

In this Appendix, we derive expressions for the bo- 
son partition and Green's functions that are essential to 
our boson and Bose-Fermi mixture AFQMC algorithms. 
These expressions have appeared in other contexts else- 
where [T2TI74) . We derive these in detail below, drawing 
from Refs. [55] and [75]. 

The fundamental relationship we aim to prove relates 
the trace of a product of one-body operators to a deter- 



minant 



e * 3 3 e » 3 3 



= Det 



(40) 



where b\ , hi are boson creation and annihilation opera- 
tors at site i and A and B are arbitrary matrices of co- 
efficients. Let us use b' to denote a row vector of boson 
creation operators: 



6 f = ••.&]*}> 



(41) 



where N is the size of the one-particle basis. Correspond- 
ingly, b will denote a column vector of annihilation oper- 
ators. A general one-body operator A is then 



i = 6t A S = ^6j A .. 6 . 



(42) 



which is a scalar and is defined by the matrix A whose 
matrix elements are given by Ajj. 



To prove Eq. 40 we first prove the following identity 

e- A e~z = e~ d , (43) 

where the matrix C defining the one- bod y operator C 

is proven, we 



-A„-B 



Once Eq. 43 



is given by e = 
can easily go to the diagonal basis to obtain Eq. [40] Let 
U'QiU = Diag[ci], where Cj are the eigenvalues of the 
matrix C, and V i — Uj^bj. Then, 



Tr h 



-KCijbj 



IIE 

mi- 



ll 



Det \\I-e 



(44) 



To prove Eq. 43 we consider the operation Aw . Using 
the boson commutation relation: bjb\ = 5jk + b\bj, we 
have 

H = £ blAijbj b\ = J2 b\A lk + blY, btA^bj, (45) 



which gives 



AV = P ■ (A + li), 



(46) 



where I is an N x N unit matrix. Note the left-hand side 
is a scalar times a row vector while the right-hand side 
is a row vector times a matrix. Repeated application of 
this equation yields 



A m & = ■ (A + iAy 



(47) 
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for any positive integer m. Thus 



= - A 6t = fet. e -(A+L4) =S t. e -A e -A 



(48) 



where in the last step the exponential can be broken up 
as the two parts commute. This is similar to the equation 
for fermions |76j . 

Now we consider an arbitrary single-boson state 



St^|o)=^^|o), 



(49) 



where <f> is a column vector containing the orbital coeffi- 
cients 4>i. The operation of the one-body propagator e~ A 
on the state leads to 

e~ A |<£) = e~ A P ■ 0|O) = V ■ e~ A ■ |0), (50) 

where in the last step we have used the fact e~ A \0) = |0). 
Similarly, for a two-boson state 

|^^) = W| ) = (gt.^)(gt.^| ) ) 

we have 

e' A (f>) = (S f • e~ A • ^)(&t • e~ A • <f>) |0). 



(51) 



(52) 



Proceeding inductively, we see that the effect of any 
single-particle propagator e~ A on any n-particle state 
(including states in which some orbitals are identical, i.e., 
multiple bosons occupying the same 1-particle orbital) is 
simply to modify each orbital by the matrix e~ A . Ap- 
plying this twice leads to the proof of Eq. |43} 

With an expression for the trace in hand, we can eval- 
uate the related boson Green's function. The Green's 
function may be written as 



G, 



Tn 


bib\e- & e- A ' 




kb]e-e~ 













(53) 



where we have used e~ A and e~ B to represent the prod- 
uct of one-boson propagators for the time slices m < k 
and m > k, respectively, with the equal-time Green's 
function measured at time-slice fc, and e~ c = e~ A e~ B . 
Transforming to the one-particle basis that diago- 

nalizes C, as in Eq. 



44 



we obtain: 



Git = 



Tr, 



Sa + 



Sij + (i 



v' 



In Tr h 



Y[e-' b >^' 



|J> 



(54) 



In equilibrium AFQMC simulations, e c represents the 
decomposition of the density matrix e~@ H as the prod- 
uct of time-sliced exponentials of quadratic operators, 
B((f>i)...B((f>i), with the corresponding time-ordering as 
defined by k, where the Green's function is measured. 
With these equations, one can readily extend fermion 
AFQMC techniques to bosons. 



APPENDIX B: ALGORITHMIC DETAILS FOR 
WORKING WITH BOSON GREEN'S 
FUNCTIONS 

The form of the boson Green's function necessitates 
three changes to the usual fermion AFQMC algorithm. 
The first two changes pertain to the equations for calcu- 
lating the ratio of determinants and the updated boson 
Green's function after each selection of a new field. The 
last pertains to the computational stability and condi- 
tioning of boson Green's functions at low temperatures. 

While the boson Green's function may be recalcu- 
lated from scratch each time it is altered, it is numer- 
ically cheaper to use the Sherman-Morrison- Woodbury 
formula, which yields the inverse of an invertible ma- 
trix plus a dyadic product. The formulas for perform- 
ing rank-one updates on the fermion Green's function 
are well-known [S31 [70]. Following Bai's derivation for 
fermions |53j , here we derive the related formulas fo r bo - 

~ c ), as given in Eq. 



54 



son Green's functions, 1/(1 — e 

Let Mi be the inverse of a boson Green's function be- 
fore the selection of a field and Mi be that after the 
selection of a field. From Eq. [54j we can write these as 

M x = I - FVi (55) 

and 



M 2 = I-FV 2 - 



(56) 



F represents a matrix appropriate for the corresponding 
C. Vi and V-z are diagonal matrices, only differing at the 
« th element. With no loss of generality, let us assume 
i = 1. Then 



where 



Vx l V 2 =I + ae 1 el, 



= VMM) 



(57) 



(58) 



As usual, e\ represents the first column of the identity 
matrix. Af 2 may then be reexpressed in terms of Mi 

M 2 = I - FVi — FVi(V^ 1 V 2 — I) 
= Mi - aFVieief 

= Mi [I + a{I - Mf^eief] . (59) 
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Expressing M 2 in terms of M\ in this form allows one 
to readily determine the ratio of determinants, r b , of the 
respective matrices. As discussed in Section III.C, r b 
must be included in the weighting factor that multiplies 
the overall walker weight after each field selection. For 
bosons, the ratio of interest is 



= Det[I/M 2 ] 
~ Det[I/Mi] 



Det[M-i] 
Det[M 2 ] 



(60) 



From above, we have 



l/r b = 



Det[M 2 ]/Det[M 1 ) 

Det[I + a{I - Mf x )eief] 

l + a(l- ef Mf 1 ei). 



Thus, 



l + a(l 



(61) 



(62) 



If one were to sample boson determinants using the 
Metropolis algorithm, it is r b that would be used in the 
acceptance criterion. 

The updated Green's function may furthermore be ob- 
tained by inverting Eq. |59| Taking the inverse, we have 



M 2 -1 = [1 + a(I - Mf^ei&f] 1 Mf l . 
Using the Sherman-Morrison- Woodbury formula, 



(A + uv 7 



A~ 



A- 1 
I 



uv T A^ 1 
v T A- 1 u 



and letting A — I , u — a(I — M 1 )ei, and v T = 
then have 



(63) 



(64) 



we 



Mf 1 



I 



m; 



a(I 



l + aef(I - M^) ei 



Mf 1 



a 



(I - Mf ^eief Mf 1 



(65) 



Since Mf 1 is simply the previous boson Green's function 
and a and r b have been calculated, this equation repre- 
sents a facile way of updating the boson Green's function. 
Analogous equations may be derived for other diagonal 
sites. 

In addition to these adjustments to the local updating 
scheme, a slight change must also be made to the way one 
inverts the boson Green's function. Just as special care 
must be taken to invert the ill-conditioned denomina- 
tor of the fermion Green's function at low temperatures, 
care must similarly be taken to invert the denominator 
of the boson Green's function. One should therefore per- 
form the same [/-DU-decomposition used for fermions [70] 
on bosons, but with a sign change reflecting the oppo- 
site sign that appears in the denominator of the boson 
Green's function: 

G b = [I -UDV}- 1 ^V-^U^V- 1 - D}- 1 ^ 1 
= V-^U DV'}- 1 !/- 1 . (66) 



In the above, U, U are orthonormal matrices, D,D' are 
diagonal matrices, and V, V are upper-triangular. 
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FIG. 1: The total, kinetic (KE), and potential (PE) energies of a 3-site Bose-Hubbard Model simulated for several values of 
Ub, tt = 0.01, and {rib} = 1 using both ED and QMC. ft denotes the inverse temperature. Agreement is within error bars for 
all points depicted. 
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FIG. 2: 3-Site Bose-Hubbard Model simulated for several values of U b , t b = 0.01, and (rib) = 1 using ED, QMC, and MFT. 
Because MFT yields the same non-interacting value of the condensate fraction regardless of Ub, only one mean-field curve is 
shown above, ft denotes the inverse temperature. Agreement between ED and QMC is exact within error bars. MFT is only 
accurate for small Ub/t. 
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FIG. 3: Number of bosons, total energies, and condensate fractions using ED, exact QMC, and the phaseless (PH) approximation 
for a 3-site Bose-Hubbard model with Ub/t = 0.5, tt = 0.01, and {rib) = 1- P denotes the inverse temperature. All points were 
produced with a timeslice of At — .025 and 50000 samples. The phaseless approximation reduces the size of the error bars on 
the number of bosons by at least half with respect to the exact QMC error bars. 
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FIG. 4: QMC vs. worm algorithm total energies for 2D Bose-Hubbard models with tb = 0.01 and (rib) = 1. /3 denotes the 
inverse temperature. Top: Total energies minus chemical potential contributions from the worm and B-AFQMC algorithms 
with decreasing temperature for a 3x3 Bose-Hubbard model for several Ub- Center: Total energies minus chemical potential 
contributions from B-AFQMC and MFT with decreasing temperature for a 3x3 Bose-Hubbard model for several Ub (note 
different scales on the horizontal axis). The QMC data is the same as used in the top panel. Bottom: Total energies with 
decreasing temperature for 2D models of varying size for Ub/t = 0.5. Total energy minus chemical potential contributions is 
plotted above in order to remove any discrepancies resulting from the fact that B-AFQMC and MFT require different chemical 
potentials to achieve the same boson densities. B-AFQMC can accurately reproduce energies for varying systems sizes and 
interaction strengths as seen by comparing to the worm algorithm. The B-AFQMC's reach is only limited by the phase problem. 
Worm and B-AFQMC energies dramatically differ from those obtained using MFT at lower temperatures. 
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FIG. 5: 2-Site Bose-Fermi-Hubbard model kinetic energies (KE), potential energies (PE), condensate fractions, and double 
occupancies per site for varying Ut — C, tt — tf = 0.01, and (nj,) = (nf) = 1 using both ED and BF-AFQMC. /3 denotes the 
inverse temperature. BF-AFQMC results are in exact agreement with those from ED. 




FIG. 6: QMC and MFT condensate fractions for the 2D Bose-Fermi-Hubbard model at U b /t = C/t = 0.5, t b = t f = 0.01, and 
(rib) = (n/) — 1. /3 denotes the inverse temperature. Good agreement is found between QMC and MFT at high temperatures. 



